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Numerical Integration of 
Stochastic Differential Equations 

By E. HELFAND 
(Manuscript received June 28, 1979) 

A procedure for numerical integration of a stochastic differential 
equation, by extension of the Runge-Kutta method, is presented. The 
technique produces results which are statistically correct to a given 
order in the time step. Second- and third-order approximations are 
explicitly displayed. 



I. INTRODUCTION 

Systematic work on numerical solution of stochastic differential 
equations (sdes) seems not to have kept pace with the considerable 
analytical developments. This parallels the lag which existed between 
the analytical and numerical study of ordinary differential equations 
near the turn of the century (which is perhaps understandable in view 
of the difficulty of implementing even straightforward algorithms at 
the time). In the last few years, there has been a burst of activity in 
performing Brownian dynamics computer simulations' to gain insight 
into motions in complex physical systems. Little attention seems to 
have been paid, though, to the systematic development of the numer- 
ical techniques in most of these works. 

In the present paper, the Runge-Kutta (rk) approximation for 
deterministic differential equations (ddes) is extended to sdes. Al- 
though we have not as yet explicitly considered other popular numer- 
ical schemes, we feel that the techniques utilized here should have 
wider applicability. For the sake of simplicity, several further restric- 
tions are placed on the discussions in this paper. These, we believe, 
can ultimately be removed by fairly simple means. 

(i) We shall work only with a single equation rather than a set of 
n equations. It has been explicitly verified that the second-order 
approximation carries over in a straightforward manner to sets (and in 
our studies of polymers 2 we used it for 600 simultaneous equations). 
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However, bear in mind Butcher's demonstration 3 that extra conditions 
arise with the rk method when one generalizes fifth-order (5o) and 
higher schemes from single equations to sets. 

(ii) We present explicit results only for low-order algorithms, sec- 
ond (2o) and third (3o) order, although the principles of higher order 
extensions will be written down. 

(Hi) Finally, we restrict attention to a simple sde, which is of the 
general form occurring in Brownian motion theory. This is 

dx/dt = f(x)+ A (t), (1) 

where the A (t) are Gaussianly distributed random variables with mean 
zero and covariance 

(A(t)A(t')) -#(*-?) (2) 

(white noise). The extension to f(x,t) appears to involve little new, but 
makes the presentation more cumbersome. 

In Section II, we review the rk technique for ddes. After defining 
more clearly what is meant by numerical solution of an sde in Sec- 
tion III, we explicitly extend the 2o RK method to sdes and outline the 
generalization to any order. A 3o RK scheme is presented in the 
appendix. Section IV is a brief discussion of the question of accuracy. 
The concluding remarks indicate areas for future studies. 

Abbreviations used throughout the paper are listed in Table I. 

II. SUMMARY OF THE RK APPROXIMATION FOR DDEs 

To set the stage, it will be useful to review 4 briefly the application of 
the rk technique to the dde 

dx/dt = f(x). (3) 

Of course, this equation can be solved by quadrature, but not when x 
and f are vectors, or when f is a function of x and t (the rk procedure 
for the latter case is presented in most standard texts 4 and does not 
differ greatly from the case we are considering). 

Begin by writing down the solution of eq. (3) as a series in the time 
step s: 

x(s) =x + sf + V4 sVo/o + (Ve)s 3 (f fo 2 + ffff )+••-, (4) 

Table I — Summary of abbreviations 

sde Stochastic differential equation 

dde Deterministic differential equation 

rk Runge-Kutta 

ko Ath order 

Is /stages 

nu, m Gaussian random variables per step 
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where f™ denotes the nth derivative of f evaluated at x„. The aim of 
many numerical procedures is to present an algorithm which, when 
expanded in s, matches the series (4) to a given order, k, in s. Of course, 
merely evaluating eq. (4) will do that, but a further aim is to avoid the 
determination of derivatives off. Thus, in the rk theory one goes from 
an initial condition x„ to x{s) in I stages by the general procedure 

g\ = f(x ), 

gi = f(Xo + /?2lS£l), 

g3 = f(x„ + fclSgi + /?32S£ 2 ), 



gi = f(x + ftnsgi + • • • + Pu-i8gi-i), (5) 

Xs = x„ + s(Aigi + A 2 g2 + • • • + Aigt). (6) 

The V2 1(1+1) parameters A\, • • • , Ai, /fei, • • • , Pu-\ are to be selected 
so that an expansion of eq. (6) in powers of s matches Eq. (4) through 
order k. Only for k < 4 can a &th order (ko) RK calculation be done in 
k stages (k s ). For k > 5, a larger number of stages than the order is 
necessary to provide enough parameters to match the true series. 
In the 2o2s RK, the parameters must satisfy 

A,+A 2 =l, (7) 

A 2 /3 21 = '/2. (8) 

This illustrates the common occurrence of less equations than param- 
eters. The user then has the freedom to select some parameters (one 
in the present case) for convenience, or to achieve the smallest error 
estimates. 5 

III. GENERALIZATION OF RK METHOD TO SDEs 

An sde does not have a definite solution. When we say that we are 
numerically integrating an sde, we mean that we are generating a 
statistically representative trajectory. Furthermore, as in numerical 
integration of a DDE, we are not going to generate the full trajectory, 
but only values of x at discrete times: x(0), x{s\), x(s\ + s 2 ), • • • . Let 
us be more specific. The stochastic process x (or set of processes) 
specified by eq. (1) is Markovian. Thus, the process is completely 
specified by the conditional probability density function p{x, s\x„), 
which gives the probability density of observing x at time s, given the 
value Xo of the variable at time zero. What we seek is a method of 
selecting a value *, with statistics correct to k\h order in s. By this, we 
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mean that the moments <x? ) are all correctly given to 0(s k )\ i.e., there 
exists a sequence C q such that for sufficiently small s 

|<*J) a -<x<s)«>e|<C,** f (9) 

for all positive integers q. The average < )„ is over the ensemble 
generated by the approximate process, while < ) e is over the exact 
process. 

An approximation algorithm will involve generation of some random 
numbers. Naturally, if p(x,s \ x a ) is known, all that need be done is to 
generate a single uniformly distributed random number, u, for each 
step and solve the equation p (x,s \x Q ) = u for x. We shall see that a lo 
approximation is equivalent to linearizing f (since f" does not enter 
until 0{s~)). For a linear sde, p(x,s\x„) is a well-known Gaussian. 6 
Use of this Gaussian as an approximate process has been suggested. 7,8 
This is practical for a single variable x, but for a large set, the amount 
of matrix manipulation is overwhelming. 

In the rk extension to be discussed, for each step of time s, 
m independent Gaussianly distributed variables, Z,(or m sets, Z,), will 
be needed. These have 

(Zi) = 0, (10) 

(ZiZj) = 8ij. (11) 

An approximation which requires m Z's will be said to be m-fold 
Gaussian, abbreviated itig. 

Now we shall present a parallel to the rk procedure for sdes. Again 
begin by developing a "power series" expansion for the solution of the 
SDE (1): 

dx/dt = f(x)+A(t). (1) 

This may be done by iteration and Taylor series expansion: 

x(s) = x + I dsifl x + I ds 2 f[x + • • • ] + w„(s\) > 

+ w„(s) (12) 

= x n + sf n + i s 2 f f + I s 3 (f f^ 2 + f 2 /JO + . . . + S, (13) 

S={w (s)) + {/Ms)} + ji/-„" J rfs,^( Sl )J 

+ j/o'Ws) + fofo[siVx(s) - W2(S)] 
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+ -f"' I dsiiv: 

b /o 



;!(*,)} 



+ ] i /yr | <fei(« - s,)ii£(s,) + i / /2" J (&,*,«£(«,) 



2' ' | o 2 



+ 5/2** J rf«i^(*i)| + •••; (14) 

w„(s)=\ dSiWn-dsi), n>0, (15) 



= r d8l <f_£^A(8 I ) f n>0 (16) 

Jo " ! 

(w„(s) is the Wiener process). The term S is a stochastic process. Its 
various parts, set off in braces, have orders in probability s 1/2 , s 3/2 , s 2 , 
s 5/2 , s 3 , • • • (N.B.: there is no s 1 term). A stochastic variable v will be 
said to have an order s k in probability if 

\(tfl)\<K q sP k , (17) 

for all positive integers q, a set of constants K q , and sufficiently 
small s. The w n are correlated Gaussian random variables with mean 
zero and covariances 

(w„(s)w m {s)) = $s n+m+1 /n\m\{n + m + 1), (18) 

(itb(8)iifa(i)) = £min(s, t). 

<*<•>*«> -ftf*-.), £2 (i9) 

The statistics of the stochastic part of the trajectory are embodied 
in the moments of S which, from eqs. (18) to (20), are 

(S) = Va s 2 £f,:' + s 3 £('/i2 flW + '/« W>" + y-u Uo v ') + ••-, (20) 

(S 2 ) = s£ + #ift + s 3 £(% ff + % UW + '/, if,',") + . . . , (21) 

<S 3 > ='/4S 3 ^" + .... (22) 

For the expansion through 0{s 3 ), the terms of S nonlinear in the w's 
do not contribute to the moments (S 4 ) and higher. Thus, these 
moments are related to the second moment by the usual Gaussian 
formulas: (S 4 ) = 3(S">, etc.; i.e., the cumulants vanish to 0(s 3 ). The 
point is that, if (S 2 ) is properly given, so will (S k ), k > 4, be. The aim 
of a ko numerical scheme will be to match not only the nonstochastic 
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terms of the series solution for x, eq. (13), but also to match all the 
moments of the stochastic term. 

We delay the presentation of the general extension of the rk 
approximation and first explicitly display a 2o2s1g scheme. Consider 
the algorithm: 

ffi-/(Xo + * 1/2 | I/2 A^), (23) 

g 2 = f(x + 8fo + s l/2 e /2 XaZ), (24) 

x = x + s (Axg x + A 2 g 2 ) + s l/2 £ 1/2 XoZ. (25) 

Z is a single Gaussian random variable with mean zero and variance 
unity, generated for each time step s. Using these equations, x can be 
developed in power of s 1/2 to 0(s 2 ): 

x = x + (Ai + A 2 )sf + A$s%W +■-■+§, (26) 

5 = \ Zs 1/2 e /2 + (AiX, + A 2 \2)Zs 3/2 Z i/2 f'o 

+ Vz(A x \\ + A^ 2 )Z 2 s\W + • • • . (27) 

The moments of 5 through O (s 2 ) are 

(§) = K(Ax\\ + AzADs^/o" + • • • (28) 

(S 2 ) = A 2 s£ + 2(A 1 \ 1 + A^XoS 2 ^ .... (29) 

To 0(s 2 ), the moments (S 3 ) and higher involve only the linear terms 
of S, so they are Gaussianly related to (S 2 ). Matching the deterministic 
part of eq. (26) to (13), eq. (28) to (20), and eq. (29) to (21), we find as 
equations for the parameters: 

Ax + A 2 = 1, (30) 

A 2 = % (31) 

\l = 1, (32) 

(AiXi + A 2 \ 2 )\„ m K, (33) 

A, A? + AzAi = x h. (34) 

The sign of Ao is immaterial since it multiplies a symmetric random 
variable. There are five equations and six parameters. A convenient 
solution set is 

A, = A 2 = x h, (35) 

- 1, (36) 

Ao = 1, (37) 

2294 THE BELL SYSTEM TECHNICAL JOURNAL, DECEMBER 1979 



and either 

(A, = 0, A 2 = 1), (38a) 

or 

(A, = 1, A 2 = 0). (38b) 

With this as background, the general procedure for constructing a 
kolsrriG approximation should be clear. Consider the m Gaussian ran- 
dom variables as a vector Z = (Z u Z 2 , • • • , Z m ). Also define Z + 1 
vectors of parameters, each of dimension m: 

A„ = (Aoi, 0, 0, ... ) 

A, = (A„,A 12 , 0, ... ) 

A/= (A n , A/2, ••• ,A, m ). (39) 

The number of scalar A parameters is m{l — Vz m + %). The general- 
ization of the rk algorithm is 

^=/(x„ + s ,/2 r /2 A,-Z), 

g* = f(x„ + sfcigi + s l/2 e /2 A 2 -Z), 

g, = f( Xo + sftug, + • • • + sfrj-igi-i + s 1/2 £ ,/2 A/-Z), (40) 

x = x„ + s(A igi + ... + A,g,) + s m l m A„-Z. (41) 

The A's and /?'s are subject to the usual rk equations since, for £ = 0, 
the dde is recovered. The equations for the A's are obtained by 
expanding eq. (41), in powers of s 1/_ ' to order s* and separating off a 
stochastic term S. Each term of the moments of § has the form of a 
product of a numerical coefficient, an integral power of s and of £, a 
product of powers of f„ and its derivatives, and a product of the A, /? 
and A parameters (the A parameters enter only as dot products of the 
A vectors). This term is equated to the term of the exact moments of 
S with the same powers of s, £, f a , and derivatives of/,, [see eqs. (20) to 
(22)]. The result is a set of equations for the A's, and the number of 
Gaussians must be chosen so that there are a sufficient number of 
parameters to satisfy these equations. One Gaussian will do for 2o2.s, 
and two Gaussians for 3o3s (see the appendix). 

IV. ACCURACY 

The accuracy of a numerical scheme for integrating a dde can be 
judged on the basis of its ability to determine trajectories for analyti- 
cally soluble equations. The schemes for sdes can only be judged on 
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a statistical basis. For example, the probability density, p(x, t),for the 
random process x defined by eq. (1) satisfies the Fokker-Planck equa- 
tion 



dp d 

dt dx 



MP-\*ip 



(42) 



This has a stationary solution 

p o (x) = N(0exp[2F(x)/Q (43) 

F(x) = I f(x')dx', (44) 



= J fixldxf, 

-I 



i/N(Q = exp \2F(x')/£\dx' (45) 

(assuming that the density is normalizable). For a stable approximation 
scheme, the distribution of x will also approach a stationary probability 
density. One could attempt to test the overall "goodness of fit" of the 
observed to the theoretical density function. 9 An easier procedure is to 
assume that eq. (43) holds and to obtain an estimate of £, for instance 
by maximum likelihood estimation. 9 The estimated £ is then compared 
with the exact £. We have used this technique and have clearly 
observed how the estimate improves with decreasing step size s. 
However, no systematic studies have been carried out yet to determine 
whether the error decreases as s k+l . 

In general, one is interested in the complete comparison of the 
transition probability p(x, s\x ) for the sde and the numerical scheme. 
This is embodied in the spectral resolution, for the exact process and 
the approximation, of p(x,s\x a ) regarded as an integral kernel. Here 
studies performed on exactly soluble systems would be of value. 

A question related to accuracy is: How long a trajectory need one 
run to reduce statistical error in some property to acceptable levels? 
The answer depends on the time, t, for decay of correlation of that 
property. New statistical information is only generated in a time of 
O(t). 10 Therefore, a simulation of total time t will lead to a decrease of 
error like (r/t) i/2 . Some systems cannot be described in such a clear- 
cut fashion since they have a spectrum of relaxation times, some of 
which may be very long. In such cases, there may be an advantage in 
reinitializing the run to break correlations. 

V. DIRECTIONS FOR FURTHER RESEARCH 

The specific procedures displayed in this paper are illustrative of 
the manner in which standard numerical techniques can be extended 
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to stochastic differential equations. There are several general direc- 
tions in which further research may be aimed. 

5.1 More general SDEs 

The numerical schemes should be directed toward more general 
sdes. The extension to sets of equations has been mentioned. More 
general forms of sdes than eq. (1) are 

dx = f(x,t)dt + 4>(x, t)dw„ (t) (46) 

or 

^ = f(x,t,A(t)). (47) 

at 

Another generalization is that A may be other than Gaussianly dis- 
tributed. Also, in the physical literature there is increased attention 
being directed to stochastic integrodifferential equations, representing 
processes with memory, such as" 12 

^ = f(x) + J drK{T)x{t - t) +A(t), (48) 



(A(t)A(t + T))ocKlr), (49) 



or more generally, 



i.i 



I 



dx 

./(*)+ | <h(;\r,xU r)\ + A, (f,0) 



with A and G related by a generalized fluctuation-dissipation theorem. 

5.2 Other numerical schemes 

It would be interesting to develop stochastic versions of other 
numerical schemes used for ddes. One may raise the objection to any 
multistep procedure that it does violence to the Markovian nature of 
the process. One would have to reuse the random variables, Z,, for 
several steps to eliminate the spurious memory to the desired order. 

5.3 General principles 

There are many matters, which are the standard fare of the deter- 
ministic numerical analyst, that should be placed in a stochastic 
context. The question of accuracy has been raised. Another is stability. 
A third question is that of step-by-step error estimation. An interesting 
problem arises in developing the analog of step-size adjustment and 
the criteria for when it is necessary. Imagine that such criteria exist 
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and a particularly large Z triggers the call for step-size adjustment. 
The new Z's that are generated should not be independent of the 
old Z's. 

Finally, as a general problem, the matter of computational speed 
should be considered. To gather statistical data, long trajectories must 
be run, sometimes on systems of many degrees of freedom. It is urgent 
that there be an analysis of various procedures with respect to their 
relative speeds, for a given accuracy. 

APPENDIX 

3 3 S 2 G Procedure 

To carry out a 3o procedure requires three stages and two Gaussian 
random variables. The explicit algorithm is eqs. (40) and (41) with 
1=3. The parameters must satisfy the equations 

A, + A2 + A3-I, (51) 

A 2 #>i + Aaifiai + #w) = % (52) 

A$i + A 8 (#,i + fa) 2 = % (53) 

Aafitafiai = %, (54) 

Aoj - 1, (55) 

Ai\n + A2A21 + AiA.,, = '/2, (56) 

A,|Xi| 8 + ;4 2 |A2|* + A 3 |A8| 2 «fc, (57) 

A,A'f, + AiAii + AjA 2 ,, = '/», (58) 

A,|A,| 2 A 11 + A. 2 |A 2 | 2 A 21 +A3|A3| 2 A3i = 1 4 (59) 

A&1A2, + A a (0 ai + B 3 2)A.„ = V* (60) 

|A,A, + A>A 2 + A 3 A 3 | 2 + 2(A 2 /32iA„ 

+ AafaiXn + i4 3 /8 32 A 2 .) = %. (61) 

The first four equations are the usual ones for a 3o RK approximation. 
They leave two degrees of freedom. A widely used solution is 

Ax = %, Ai = %, A = %; (62) 

/? 21 = '/2, ftl = 0, /?32 = %. (63) 

With this set, the remaining seven equations can be solved for the 
A parameters. The solution is 
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X.i-1, A„ = 0, \ 2 i = ^ A 3 , = %. (64) 

There are four solutions for the X/ 2s , two of which are complex. The 
real solutions are either 

\i2 - 0.245538, 
\ 22 = -0.023225, 

X„ = 0.544169, (65) 

or \i2 = -1.34583, 

\ 22 = 1.24987, 
X.,2 = 0.385032. (66) 

Solution (65) is probably superior because it uses less of the Z 2 process. 
(All the An and/or all the A /2 may be reversed in sign as an acceptable 
solution, as is evident since they multiply symmetrically distributed 
random numbers.) 
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